Ground state properties of a homogeneous 
2D system of Bosons with dipolar interactions. 
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The ground-state phase properties of a two-dimensional Bose system with dipole-dipole interac- 
tions is studied by means of quantum Monte Carlo techniques. Limitations of mean-field theory in 
a two-dimensional geometry are discussed. A quantum phase transition from gas to solid is found. 
Crystal is tested for existence of a supersolid in the vicinity of the phase transition. Existence of 
mesoscopic analogue of the off-diagonal long-range order is shown in the one-body density matrix 
in a finite-size crystal. Non-zero superfluid fraction is found in a finite-size crystal, the signal being 
dramatically increased in presence of vacancies. 
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I. MODEL AND METHODS 

We study properties of a two-dimensional (2D) system of bosons with dipolar interaction. We consider a polarized 
system and assume that dipolar moments are oriented perpendicularly to the 2D plane. This assures that the 



interaction potential Vint (f ) = C, 



ddl 



is always repulsive and there are no instabilities caused by dipolar attraction. 



The following model Hamiltonian is used to describe the system: 
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where m is the dipolar mass and N the number of dipoles. 

Properties of a homogeneous system are governed by one characteristic parameter, the dimensionless density nrp, 
where the characteristic length rg is proportional to the interaction strength: tq = mCdd/ 4:7711^ . Deeply in the 
dilute regime one expects that a short-range interaction potential (like the dipolar one) can be described by only one 
parameter, namely, the s-wave scattering length a. Parameters rp and a are directly related: a = e^'^rg = 3.17...ro. 

We perform a numerical study of the ground-state properties of this system using quantum Monte Carlo methods. 
Firstly, using the variational Monte Carlo (VMC) method it is possible to evaluate multidimensional averages over 
the trial wavefunction ipx- In the calculation of the energy and superfluid fraction ns/n, the variational parameters 
in ^jt are chosen such that they minimize variational energy. In the calculation of the one-body density matrix we 
optimize parameters so that the difference between variational and mixed estimators is minimal. Secondly, we use the 
diffusion Monte Carlo (DMC) method based on solving Schrodinger equation in imaginary time at T = 0. The DMC 
method permits to find the ground state energy i? of a bosonic system exactly (in statistical sense) . Also superfluid 




FIG. 1: Correlation functions in gas phase: Pair distribution function. 



FIG. 2: Correlation functions in gas phase: One-body density matrix. 



density ns/n and local quantities (e.g. 52 (-z), Sk, etc.) can be found in a "pure" (non-depending on the choice of trial 
w.f.) way. An extrapolation procedure can be used for predictions of non-local quantities (e.g. gi{z), rik, etc.). 

II. TRIAL WAVE FUNCTION 

We construct the trial wave function (w.f.) in the following form: 

N M / N \ 

^^(r-i, r^^) = J] f2{\n " ^.1) x [] E exp{-«(n - r^O^ (2) 

i<j k=l I 

where r^jZ = 1,A^ are particle coordinates and r^^^^- ,k — l,Af are coordinates of triangular lattice sites. Trial w.f. 
^ is symmetric under exchange of any two particles. The two-body Jastrow term ^i(r) is chosen[ll 0] at short 
distances as a solution of the 2-body scattering problem at zero energy. At large distances, the functional form of the 
hydrodynamic solution is used Thus, /2(?') accounts for pair-collisions relevant for short distances and collective 
behavior (sound) at large distances. One should note that V't is not of a Nosanow-Jastrow type, as by moving one 
particle all M terms in the product in are changed, thus introducing a global change (i.e. depending as well on 
coordinates of other particles), so that this term is not a one-body term, but rather a many-body term. Another 
feature of this w.f. is that number of the particles TV can be different from number of lattice sites M making it suitable 
for studying of a crystal with vacancies. 

The parameter a describes particle localization close to lattice sites. Typical dependence (at sufficiently large 
density, nr\ ^ 10) of variational energy on parameter a shows two minima. Position of the first minimum is a = 0. 
In this case, the translational invariance is preserved and density profile is flat. This minimum corresponds to a 
gas/liquid state. In the second minimum a is finite and translational invariance is broken. Density profile has crystal 
symmetry. This minimum corresponds to a solid state. Thus, with the single trial w.f. and different variational 
parameters we are able to describe distinct phases. 

III. RESULTS 
A. Quantum phase transition 

The ground state phase at small densities corresponds to a gas, solid being metastable. On the contrary, at large 
densities solid state is energetically preferable. The critical density of the quantum phase transition ncr\ — 290(30) 
was obtained by constructing fits to the energy of gas and solid phases ^2] . This estimation of the critical density 
is in agreement with a Path Integral Monte Carlo calculation [1] done at low finite temperature Ucr^ — 320(140). 
Green's Function Monte Carlo calculation [5] provided a slightly lower critical density UfT^ = 230(20) but in this case 
a discrete model was used, and by decreasing the filling factor, a small increase in the critical value was found. 

Fig. [H shows the pair distribution function g2(j) — (5'^(0)^'t(r)^'(r)^'(0))/n^ in the gas phase in a wide range of 
densities. At the largest density (close to the phase transition) there is a well pronounced first peak followed by a 
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FIG. 3: Solid phase, nrg = 290: one-body density matrix for different system sizes. Solid line — non symmetrized w.f. 



number of well visible oscillations. This is a manifestation of strong correlations present in the system close to the 
point of a quantum phase transition. As the density is lowered the height of the peak is decreased and eventually it 
disappears, leading to a smooth behavior without any visible oscillations. This smooth behavior is characteristic for 
weakly-interacting Bose systems. 

B. Dilute regime 

In the dilute regime one expects the mean-field theory to be applicable. As was rigorously derived in Ref. [ 6]], the 
2D Gross-Pitaevskii equation (GPE) has a coupling constant g2D — Aufi^ /m\ ln(na^)| dependent on the density. This 
leads to a logarithmic dependence of the mean-field energy on the density [3] 

Emf 1 27r?i2 

It turns out that the mean-field contribution to the energy ^ is the only well established term. We perform a study 
of beyond mean-field terms. Comparison to numerical results for hard-disks[8] show that for densities na^ < 10~^ 
the exact shape of the interaction potential is no longer important, and that the only relevant parameter is s-wave 
scattering length. The peculiarity of a two-dimensional system is that for such a small densities, 10~^° < na^ < 10~^, 
there is a notable difference (of several percent) between MF-GPE and exact result. Moreover, to our knowledge 
there is no analytical theory, able to reproduce correctly the energy in the whole region of the universal regime. 
Summarizing, the mean-field description has limitations (failure) in a two-dimensional system in the universal regime. 

The one-body density matrix gi{r) — (^'T(r)^(0))/n in the gas phase has a finite asymptotic value, as reported in 
Fig. [21 In the thermodynamic limit (TV —* oo) finite asymptotic values are manifestations of off-diagonal long-range 
order (ODLRO). The asymptotic value gives the condensate fraction. In the dilute regime almost all of the particles 
are condensed, but increasing the density the stronger interactions deplete the condensate and condensation fraction 
drops down to 1% close to the phase transition point. An important question is what happens to the condensate as 
system crystallizes. 

C. Study of a supersolid 

There are several definitions of a supersolid: 

1. Spatial order of a solid + finite superfluid density 

2. Spatial order of a solid (broken- symmetry oscillations in diagonal element of OBDM) + off-diagonal long-range 
order in OBDM 

Generally, it is believed that both definitions are equivalent. 
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FIG. 4: Solid phase, nrg — 290: winding number, open symbols: N — M = 12 (circles), 16 (diamonds), 24 (squares); solid 
symbols M = 24,7V — 23 (lower curve), M = 24; = 22 (upper curve); dashed line — non symmetrized w.f; dotted line — 
diffusion constant of a free particle. 



Reduced dimensionality increases the role of quantum fluctuations. This makes a two-dimensional crystal a good 
candidate for having a supersolid. Notice that in a one-dimensional, system quantum fluctuations destroy crystalline 
long-range diagonal order. A previous study for the presence a supersolid in a two-dimensional dipolar system is not 
conclusive. Low-temperature (PIMC) simulation^] shows that gas phase is completely superfluid, while no superfluid 
fraction is found in crystal phase. Still, the presence of (a possible) supersolid can be masked by much smaller critical 
temperature in a crystal. A zero-temperature method was used with a symmetrized trial w.f. in [[^]. No conclusions 
were drawn for presence/absence of a supersolid due to an unsufficient overlap of the trial w.f. with the actual ground 
state. 

Fig. [3] shows the one-body density matrix gi{r) in the crystal phase close to the phase transition. While energy 
in a crystal is not sensitive to symmetrization of the wavefunction, it is crucial to symmetrize w.f. in the calculation 
of gi{r). Indeed, without symmetrization off-diagonal element gi{r) decays exponentially fast to zero, see thick line 
in Fig. [31 Using symmetrized w.f. we find instead a finite asymptotic value (of the order of 3 x 10"'' for N = 108 
particles). There is a certain decay of the condensate fraction as system size increases. 

Finite-size effects are very important in a 2D dipolar system. Indeed, the characteristic dependence of the energy, 
OBDM limiting value, etc. is 1/\/N instead of 1/iV in short-range potentials as can be seen from the tail correction 
of the potential energy. For this reason a proper study of the supersolid in the thermodynamic limit should be done. 
Here, we limit ourselves to some preliminary results for a finite-size system. 

The superfluid fraction corresponds to the slope of the winding number (diffusion coefficient D of the center of 
masses in imaginary time r). Once again, symmetrization is crucial, otherwise artificial zero slope is obtained (dashed 
line in Fig. |4]). Using properly symmetrized w.f. we find superfluid signal. The signal gets weaker as we increase 
number of particles N (see open symbols in Fig. U]). 

The trial w.f. ^ permits us to investigate role of vacancies. A system with M = 24 lattice sites and 0; 1; 2 vacancies 
is studied. The superfluid density experiences dramatic effects in the presence of vacancies. The signal increases from 
w 2% for vacancies to « 40% for 2 vacancies. 



IV. CONCLUSIONS 

To summarize, the Diffusion Monte Carlo method was used to study the properties of a dipolar 2D Bose system 
at r = 0. The ground state energy, pair distribution function, one-body density matrix were calculated in a wide 
range of densities. The gas-solid quantum phase transition is found at density nr^ = 290(30). Limitations (failure) of 
mean-field description were pointed out in the universal low-density regime. Existence of mesoscopic analogue of the 
off-diagonal long-range order was shown in one-body density matrix in a finite-size crystal close to phase transition. 
Non-zero superfluid fraction was found in a finite-size crystal. The superfluid signal is dramatically increased in 
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presence of vacancies. 
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